Compiled under R version 4.2.0 (2022-04-22)

WARNING: edit the working directory to your preferred folder.

This document details all analyses performed in R for the study:
Legendre, L. J., S. Choi, and J. A. Clarke. The use of inconsistent terminology for reptile eggshell traits affects the outcome of evolutionary analyses. Journal of Anatomy.

For more information regarding the study, datasets, and analyses, please refer to the Main Text and Supplementary Information of this paper. If you have any additional questions, feel free to email me at .

Loading packages

library(ape)
library(castor)
library(evobiR)
library(ggtree)
library(phytools)
library(RColorBrewer)

Data and tree

  • Tree
tree<-read.nexus("treewhole_newversion.trees.nex")
plotTree(tree, fsize=0.5,lwd=1,type="fan",ftype="i")

  • Data
data<-read.table("Datawhole_newproject.txt", header=T)
setdiff(tree$tip.label, data$Taxon) # taxa and data match
## character(0)
rownames(data)<-data$Taxon
datanew<-ReorderData(tree, data)

Ancestral state reconstruction: discrete trait (soft/semi-rigid/hard)

  • Extract data vectors for each coding
thisstudy<-datanew$Type_2021; names(thisstudy)<-rownames(datanew)
norell<-datanew$Type_Norell; names(norell)<-rownames(datanew)
legendre<-datanew$Type_Legendre; names(legendre)<-rownames(datanew)

# Color palette
cols<-setNames(c("royalblue","green3","red3"),c("Soft","Semi-rigid","Hard"))

# Visualize tree with node numbers
ggtree(tree) + geom_text2(aes(subset=!isTip, label=node), hjust=-.3) + geom_tiplab()

Using character coding defined in this study (new scoring)

  • Perform SIMMAP using phytools – 1000 iterations, using AIC to select best model (code modified from Liam Revell, see here)
x<-thisstudy

aic<-function(logL,k) 2*k-2*logL
aic.w<-function(aic){
    d.aic<-aic-min(aic)
    exp(-1/2*d.aic)/sum(exp(-1/2*d.aic))
}

logL<-sapply(c("ER","SYM","ARD"),
    function(model,tree,x) make.simmap(tree,x,model)$logL,
    tree=tree,x=x)
## make.simmap is sampling character histories conditioned on
## the transition matrix
## 
## Q =
##                    Hard   Semi-rigid         Soft
## Hard       -0.001035456  0.000517728  0.000517728
## Semi-rigid  0.000517728 -0.001035456  0.000517728
## Soft        0.000517728  0.000517728 -0.001035456
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##       Hard Semi-rigid       Soft 
##  0.3333333  0.3333333  0.3333333
## Done.
## make.simmap is sampling character histories conditioned on
## the transition matrix
## 
## Q =
##                     Hard    Semi-rigid          Soft
## Hard       -0.0010899924  0.0004565271  0.0006334653
## Semi-rigid  0.0004565271 -0.0007161171  0.0002595901
## Soft        0.0006334653  0.0002595901 -0.0008930554
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##       Hard Semi-rigid       Soft 
##  0.3333333  0.3333333  0.3333333
## Done.
## make.simmap is sampling character histories conditioned on
## the transition matrix
## 
## Q =
##                    Hard    Semi-rigid          Soft
## Hard       0.0000000000  0.0000000000  0.0000000000
## Semi-rigid 0.0072145631 -0.0109014209  0.0036868578
## Soft       0.0004663626  0.0002355611 -0.0007019237
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##       Hard Semi-rigid       Soft 
##  0.3333333  0.3333333  0.3333333
## Done.
logL
##        ER       SYM       ARD 
## -56.68015 -56.39465 -45.61188
AIC<-mapply(aic,logL,c(1,3,6))
AIC
##       ER      SYM      ARD 
## 115.3603 118.7893 103.2238
AIC.W<-aic.w(AIC)
AIC.W
##           ER          SYM          ARD 
## 0.0023088700 0.0004157184 0.9972754116
nsim<-1000
Nsim<-round(nsim*AIC.W)
d<-if(sum(Nsim)>nsim) -1 else 1
nsim<-Nsim+d*sample(c(rep(1,abs(nsim-sum(Nsim))),
    rep(0,length(Nsim)-abs(nsim-sum(Nsim)))))
nsim
##  ER SYM ARD 
##   2   1 997
o1<-make.simmap(tree,x,model="ARD",nsim=nsim[1])
## make.simmap is sampling character histories conditioned on
## the transition matrix
## 
## Q =
##                    Hard    Semi-rigid          Soft
## Hard       0.0000000000  0.0000000000  0.0000000000
## Semi-rigid 0.0072145631 -0.0109014209  0.0036868578
## Soft       0.0004663626  0.0002355611 -0.0007019237
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##       Hard Semi-rigid       Soft 
##  0.3333333  0.3333333  0.3333333
## Done.
o2<-make.simmap(tree,x,model="ER",nsim=nsim[2])
## make.simmap is sampling character histories conditioned on
## the transition matrix
## 
## Q =
##                    Hard   Semi-rigid         Soft
## Hard       -0.001035456  0.000517728  0.000517728
## Semi-rigid  0.000517728 -0.001035456  0.000517728
## Soft        0.000517728  0.000517728 -0.001035456
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##       Hard Semi-rigid       Soft 
##  0.3333333  0.3333333  0.3333333
## Done.
treesstudy<-c(o1,o2)

objstudy<-describe.simmap(treesstudy)

plot(objstudy,type="fan",fsize=0.01,lwd=1,ftype="i", colors=cols,ylim=c(-2,Ntip(tree)),offset=20, part=0.97)
add.simmap.legend(colors=cols,prompt=FALSE,x=-300,y=280)

Here, a semi-rigid eggshell is the ancestral state for almost all major reptile clades – reptiles, lepidosaurs, squamates, turtles, archelosaurs, archosaurs, dinosaurs,ornithischians, and saurischians. However, there are only 6 taxa with semi-rigid eggshells out of 208!

=> Some of these taxa may have a excessive influence on this result – probably the two sauropodomorphs with that eggshell type, since they are the closer in age to the root and to the aforementioned subclades.

Using character coding from Norell et al. (2020) – ratio scoring

x<-norell

logL<-sapply(c("ER","SYM","ARD"),
    function(model,tree,x) make.simmap(tree,x,model)$logL,
    tree=tree,x=x)
## make.simmap is sampling character histories conditioned on
## the transition matrix
## 
## Q =
##                     Hard    Semi-rigid          Soft
## Hard       -0.0013917968  0.0006958984  0.0006958984
## Semi-rigid  0.0006958984 -0.0013917968  0.0006958984
## Soft        0.0006958984  0.0006958984 -0.0013917968
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##       Hard Semi-rigid       Soft 
##  0.3333333  0.3333333  0.3333333
## Done.
## make.simmap is sampling character histories conditioned on
## the transition matrix
## 
## Q =
##                     Hard    Semi-rigid          Soft
## Hard       -0.0013792703  0.0004831551  0.0008961152
## Semi-rigid  0.0004831551 -0.0010895186  0.0006063635
## Soft        0.0008961152  0.0006063635 -0.0015024787
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##       Hard Semi-rigid       Soft 
##  0.3333333  0.3333333  0.3333333
## Done.
## make.simmap is sampling character histories conditioned on
## the transition matrix
## 
## Q =
##                   Hard   Semi-rigid         Soft
## Hard       0.000000000  0.000000000  0.000000000
## Semi-rigid 0.008810469 -0.014180320  0.005369851
## Soft       0.001156338  0.001443489 -0.002599828
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##       Hard Semi-rigid       Soft 
##  0.3333333  0.3333333  0.3333333
## Done.
logL
##        ER       SYM       ARD 
## -70.17246 -69.66276 -57.83996
AIC<-mapply(aic,logL,c(1,3,6))
AIC
##       ER      SYM      ARD 
## 142.3449 145.3255 127.6799
AIC.W<-aic.w(AIC)
AIC.W
##           ER          SYM          ARD 
## 0.0006534088 0.0001472160 0.9991993752
nsim<-1000
Nsim<-round(nsim*AIC.W)
d<-if(sum(Nsim)>nsim) -1 else 1
nsim<-Nsim+d*sample(c(rep(1,abs(nsim-sum(Nsim))),
    rep(0,length(Nsim)-abs(nsim-sum(Nsim)))))
nsim
##  ER SYM ARD 
##   1   0 999
n1<-make.simmap(tree,x,model="ARD",nsim=nsim[1])
## make.simmap is sampling character histories conditioned on
## the transition matrix
## 
## Q =
##                   Hard   Semi-rigid         Soft
## Hard       0.000000000  0.000000000  0.000000000
## Semi-rigid 0.008810469 -0.014180320  0.005369851
## Soft       0.001156338  0.001443489 -0.002599828
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##       Hard Semi-rigid       Soft 
##  0.3333333  0.3333333  0.3333333
## Done.
n2<-make.simmap(tree,x,model="ER",nsim=nsim[2])
## make.simmap is sampling character histories conditioned on
## the transition matrix
## 
## Q =
##                     Hard    Semi-rigid          Soft
## Hard       -0.0013917968  0.0006958984  0.0006958984
## Semi-rigid  0.0006958984 -0.0013917968  0.0006958984
## Soft        0.0006958984  0.0006958984 -0.0013917968
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##       Hard Semi-rigid       Soft 
##  0.3333333  0.3333333  0.3333333
## Done.
treesnorell<-c(o1,o2)

objnorell<-describe.simmap(treesnorell)

plot(objnorell,type="fan",fsize=0.01,lwd=1,ftype="i",
     colors=cols,ylim=c(-2,Ntip(tree)),part=0.97)
add.simmap.legend(colors=cols,prompt=FALSE,x=-300,y=280)

Completely different result – let us check which taxa are different from the previous ASR…

datanew[c(which(datanew$Type_2021!=datanew$Type_Norell)),]
# Only 7 taxa are different – two non-avian dinosaurs and five turtles

# Changing character state for the two non-avian dinosaurs
x[c(21,22)]<-"Semi-rigid"

logL<-sapply(c("ER","SYM","ARD"),
    function(model,tree,x) make.simmap(tree,x,model)$logL,
    tree=tree,x=x)
## make.simmap is sampling character histories conditioned on
## the transition matrix
## 
## Q =
##                     Hard    Semi-rigid          Soft
## Hard       -0.0014141228  0.0007070614  0.0007070614
## Semi-rigid  0.0007070614 -0.0014141228  0.0007070614
## Soft        0.0007070614  0.0007070614 -0.0014141228
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##       Hard Semi-rigid       Soft 
##  0.3333333  0.3333333  0.3333333
## Done.
## make.simmap is sampling character histories conditioned on
## the transition matrix
## 
## Q =
##                     Hard    Semi-rigid          Soft
## Hard       -0.0013866412  0.0005855999  0.0008010413
## Semi-rigid  0.0005855999 -0.0013144727  0.0007288728
## Soft        0.0008010413  0.0007288728 -0.0015299141
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##       Hard Semi-rigid       Soft 
##  0.3333333  0.3333333  0.3333333
## Done.
## make.simmap is sampling character histories conditioned on
## the transition matrix
## 
## Q =
##                    Hard    Semi-rigid         Soft
## Hard       0.0000000000  0.0000000000  0.000000000
## Semi-rigid 0.0081716323 -0.0131662651  0.004994633
## Soft       0.0005407866  0.0007809644 -0.001321751
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##       Hard Semi-rigid       Soft 
##  0.3333333  0.3333333  0.3333333
## Done.
AIC<-mapply(aic,logL,c(1,3,6))
AIC.W<-aic.w(AIC)
nsim<-1000
Nsim<-round(nsim*AIC.W)
d<-if(sum(Nsim)>nsim) -1 else 1
nsim<-Nsim+d*sample(c(rep(1,abs(nsim-sum(Nsim))),
    rep(0,length(Nsim)-abs(nsim-sum(Nsim)))))
nsim
##   ER  SYM  ARD 
##    0    0 1000
treesnorell<-make.simmap(tree,x,model="ARD",nsim=1000)
## make.simmap is sampling character histories conditioned on
## the transition matrix
## 
## Q =
##                    Hard    Semi-rigid         Soft
## Hard       0.0000000000  0.0000000000  0.000000000
## Semi-rigid 0.0081716323 -0.0131662651  0.004994633
## Soft       0.0005407866  0.0007809644 -0.001321751
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##       Hard Semi-rigid       Soft 
##  0.3333333  0.3333333  0.3333333
## Done.
objnorell<-describe.simmap(treesnorell)
plot(objnorell,fsize=0.5,lwd=1,ftype="i",colors=cols,ylim=c(-2,Ntip(tree)),part=0.95)
add.simmap.legend(colors=cols,prompt=FALSE,x=20,y=180)

All major clades are now recovered as semi-rigid again => strong influence of the two sauropodomorphs.

Using character coding from Legendre et al. (2020) – shell unit scoring

x<-legendre

logL<-sapply(c("ER","SYM","ARD"),
    function(model,tree,x) make.simmap(tree,x,model)$logL,
    tree=tree,x=x)
## make.simmap is sampling character histories conditioned on
## the transition matrix
## 
## Q =
##                     Hard    Semi-rigid          Soft
## Hard       -0.0007487292  0.0003743646  0.0003743646
## Semi-rigid  0.0003743646 -0.0007487292  0.0003743646
## Soft        0.0003743646  0.0003743646 -0.0007487292
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##       Hard Semi-rigid       Soft 
##  0.3333333  0.3333333  0.3333333
## Done.
## make.simmap is sampling character histories conditioned on
## the transition matrix
## 
## Q =
##                     Hard    Semi-rigid          Soft
## Hard       -0.0006125127  0.0000000000  0.0006125127
## Semi-rigid  0.0000000000 -0.0005646183  0.0005646183
## Soft        0.0006125127  0.0005646183 -0.0011771310
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##       Hard Semi-rigid       Soft 
##  0.3333333  0.3333333  0.3333333
## Done.
## make.simmap is sampling character histories conditioned on
## the transition matrix
## 
## Q =
##                     Hard   Semi-rigid          Soft
## Hard       -0.0004334394  0.000000000  0.0004334394
## Semi-rigid  0.0034884449 -0.007010172  0.0035217272
## Soft        0.0003883791  0.000000000 -0.0003883791
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##       Hard Semi-rigid       Soft 
##  0.3333333  0.3333333  0.3333333
## Done.
logL
##        ER       SYM       ARD 
## -45.68062 -43.28735 -41.37627
AIC<-mapply(aic,logL,c(1,3,6))
AIC
##       ER      SYM      ARD 
## 93.36124 92.57469 94.75254
AIC.W<-aic.w(AIC)
AIC.W
##        ER       SYM       ARD 
## 0.3355056 0.4971605 0.1673339
nsim<-1000
Nsim<-round(nsim*AIC.W)
d<-if(sum(Nsim)>nsim) -1 else 1
nsim<-Nsim+d*sample(c(rep(1,abs(nsim-sum(Nsim))),
    rep(0,length(Nsim)-abs(nsim-sum(Nsim)))))
nsim
##  ER SYM ARD 
## 336 497 167
l1<-make.simmap(tree,x,model="ER",nsim=nsim[1])
## make.simmap is sampling character histories conditioned on
## the transition matrix
## 
## Q =
##                     Hard    Semi-rigid          Soft
## Hard       -0.0007487292  0.0003743646  0.0003743646
## Semi-rigid  0.0003743646 -0.0007487292  0.0003743646
## Soft        0.0003743646  0.0003743646 -0.0007487292
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##       Hard Semi-rigid       Soft 
##  0.3333333  0.3333333  0.3333333
## Done.
l2<-make.simmap(tree,x,model="SYM",nsim=nsim[2])
## make.simmap is sampling character histories conditioned on
## the transition matrix
## 
## Q =
##                     Hard    Semi-rigid          Soft
## Hard       -0.0006125127  0.0000000000  0.0006125127
## Semi-rigid  0.0000000000 -0.0005646183  0.0005646183
## Soft        0.0006125127  0.0005646183 -0.0011771310
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##       Hard Semi-rigid       Soft 
##  0.3333333  0.3333333  0.3333333
## Done.
l3<-make.simmap(tree,x,model="ARD",nsim=nsim[3])
## make.simmap is sampling character histories conditioned on
## the transition matrix
## 
## Q =
##                     Hard   Semi-rigid          Soft
## Hard       -0.0004334394  0.000000000  0.0004334394
## Semi-rigid  0.0034884449 -0.007010172  0.0035217272
## Soft        0.0003883791  0.000000000 -0.0003883791
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##       Hard Semi-rigid       Soft 
##  0.3333333  0.3333333  0.3333333
## Done.
treeslegendre<-c(l1,l2,l3)

objlegendre<-describe.simmap(treeslegendre)

plot(objlegendre,type="fan",fsize=0.01,lwd=1,ftype="i",
     colors=cols,ylim=c(-2,Ntip(tree)),part=0.97)
add.simmap.legend(colors=cols,prompt=FALSE,x=-300,y=280)

Reptiles, lepidosaurs, and squamates are recovered as ancestrally soft-shelled, while turtles, archelosaurs, archosaurs, and all dinosaur clades are recovered as hard-shelled.

The pattern seems to follow whichever character state is coded for the two sauropodomorphs Lufengosaurus and Massospondylus.

To test this hypothesis, we must remove all other taxa (n = 7) that are not coded the same in all three studies, and see if the pattern is still present.

  • Remove all taxa with differences in coding, except the two sauropodomorphs
remove<-rownames(datanew[c(which(datanew$Type_Legendre!=datanew$Type_Norell)),][c(3:9),])
x<-thisstudy[!names(thisstudy)%in%remove]
treenew<-drop.tip(tree, setdiff(tree$tip.label, names(x)))

# Visualize tree with node numbers
ggtree(treenew) + geom_text2(aes(subset=!isTip, label=node), hjust=-.3) + geom_tiplab()

Replicate analysis on the new tree 3 times, with identical coding for all taxa except the two sauropodomorphs

  • Coded as semi-rigid
logL<-sapply(c("ER","SYM","ARD"),
    function(model,treenew,x) make.simmap(treenew,x,model)$logL,
    tree=treenew,x=x)
## make.simmap is sampling character histories conditioned on
## the transition matrix
## 
## Q =
##                     Hard    Semi-rigid          Soft
## Hard       -0.0008797485  0.0004398742  0.0004398742
## Semi-rigid  0.0004398742 -0.0008797485  0.0004398742
## Soft        0.0004398742  0.0004398742 -0.0008797485
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##       Hard Semi-rigid       Soft 
##  0.3333333  0.3333333  0.3333333
## Done.
## make.simmap is sampling character histories conditioned on
## the transition matrix
## 
## Q =
##                     Hard    Semi-rigid          Soft
## Hard       -7.696326e-04  6.160148e-05  0.0007080311
## Semi-rigid  6.160148e-05 -7.755548e-04  0.0007139533
## Soft        7.080311e-04  7.139533e-04 -0.0014219844
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##       Hard Semi-rigid       Soft 
##  0.3333333  0.3333333  0.3333333
## Done.
## make.simmap is sampling character histories conditioned on
## the transition matrix
## 
## Q =
##                     Hard    Semi-rigid          Soft
## Hard       -0.0002017335  0.0000000000  0.0002017335
## Semi-rigid  0.0037018375 -0.0037018375  0.0000000000
## Soft        0.0020196421  0.0008852088 -0.0029048508
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##       Hard Semi-rigid       Soft 
##  0.3333333  0.3333333  0.3333333
## Done.
AIC<-mapply(aic,logL,c(1,3,6))
AIC.W<-aic.w(AIC)

nsim<-1000
Nsim<-round(nsim*AIC.W)
d<-if(sum(Nsim)>nsim) -1 else 1
nsim<-Nsim+d*sample(c(rep(1,abs(nsim-sum(Nsim))),
    rep(0,length(Nsim)-abs(nsim-sum(Nsim)))))
nsim
##  ER SYM ARD 
## 388 195 417
s1<-make.simmap(treenew,x,model="ER",nsim=nsim[1])
## make.simmap is sampling character histories conditioned on
## the transition matrix
## 
## Q =
##                     Hard    Semi-rigid          Soft
## Hard       -0.0008797485  0.0004398742  0.0004398742
## Semi-rigid  0.0004398742 -0.0008797485  0.0004398742
## Soft        0.0004398742  0.0004398742 -0.0008797485
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##       Hard Semi-rigid       Soft 
##  0.3333333  0.3333333  0.3333333
## Done.
s2<-make.simmap(treenew,x,model="SYM",nsim=nsim[2])
## make.simmap is sampling character histories conditioned on
## the transition matrix
## 
## Q =
##                     Hard    Semi-rigid          Soft
## Hard       -7.696326e-04  6.160148e-05  0.0007080311
## Semi-rigid  6.160148e-05 -7.755548e-04  0.0007139533
## Soft        7.080311e-04  7.139533e-04 -0.0014219844
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##       Hard Semi-rigid       Soft 
##  0.3333333  0.3333333  0.3333333
## Done.
s3<-make.simmap(treenew,x,model="ARD",nsim=nsim[3])
## make.simmap is sampling character histories conditioned on
## the transition matrix
## 
## Q =
##                     Hard    Semi-rigid          Soft
## Hard       -0.0002017335  0.0000000000  0.0002017335
## Semi-rigid  0.0037018375 -0.0037018375  0.0000000000
## Soft        0.0020196421  0.0008852088 -0.0029048508
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##       Hard Semi-rigid       Soft 
##  0.3333333  0.3333333  0.3333333
## Done.
treenewsr<-c(s1,s2,s3)

objsr<-describe.simmap(treenewsr)

plot(objsr,type="fan",fsize=0.01,lwd=1,ftype="i",
     colors=cols,ylim=c(-2,Ntip(treenew)),part=0.97)
add.simmap.legend(colors=cols,prompt=FALSE,x=-300,y=280)

Most inner clades are now recovered as soft-shelled – most likely an artefact due to the position of Mussaurus – closer to the other inner nodes than any other taxon, since it is the oldest specimen in the sample.

  • Coded as soft
x[c(21,22)]<-"Soft"

logL<-sapply(c("ER","SYM","ARD"),
    function(model,treenew,x) make.simmap(treenew,x,model)$logL,
    tree=treenew,x=x)
## make.simmap is sampling character histories conditioned on
## the transition matrix
## 
## Q =
##                     Hard    Semi-rigid          Soft
## Hard       -0.0008695627  0.0004347813  0.0004347813
## Semi-rigid  0.0004347813 -0.0008695627  0.0004347813
## Soft        0.0004347813  0.0004347813 -0.0008695627
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##       Hard Semi-rigid       Soft 
##  0.3333333  0.3333333  0.3333333
## Done.
## make.simmap is sampling character histories conditioned on
## the transition matrix
## 
## Q =
##                     Hard    Semi-rigid          Soft
## Hard       -0.0007415283  0.0000000000  0.0007415283
## Semi-rigid  0.0000000000 -0.0005402947  0.0005402947
## Soft        0.0007415283  0.0005402947 -0.0012818231
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##       Hard Semi-rigid       Soft 
##  0.3333333  0.3333333  0.3333333
## Done.
## make.simmap is sampling character histories conditioned on
## the transition matrix
## 
## Q =
##                     Hard    Semi-rigid          Soft
## Hard       -0.0001997071  0.0000000000  0.0001997071
## Semi-rigid  0.0023041591 -0.0023041591  0.0000000000
## Soft        0.0021232656  0.0005497209 -0.0026729866
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##       Hard Semi-rigid       Soft 
##  0.3333333  0.3333333  0.3333333
## Done.
AIC<-mapply(aic,logL,c(1,3,6))
AIC.W<-aic.w(AIC)

nsim<-1000
Nsim<-round(nsim*AIC.W)
d<-if(sum(Nsim)>nsim) -1 else 1
nsim<-Nsim+d*sample(c(rep(1,abs(nsim-sum(Nsim))),
    rep(0,length(Nsim)-abs(nsim-sum(Nsim)))))
nsim
##  ER SYM ARD 
## 113 284 603
so1<-make.simmap(treenew,x,model="ER",nsim=nsim[1])
## make.simmap is sampling character histories conditioned on
## the transition matrix
## 
## Q =
##                     Hard    Semi-rigid          Soft
## Hard       -0.0008695627  0.0004347813  0.0004347813
## Semi-rigid  0.0004347813 -0.0008695627  0.0004347813
## Soft        0.0004347813  0.0004347813 -0.0008695627
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##       Hard Semi-rigid       Soft 
##  0.3333333  0.3333333  0.3333333
## Done.
so2<-make.simmap(treenew,x,model="SYM",nsim=nsim[2])
## make.simmap is sampling character histories conditioned on
## the transition matrix
## 
## Q =
##                     Hard    Semi-rigid          Soft
## Hard       -0.0007415283  0.0000000000  0.0007415283
## Semi-rigid  0.0000000000 -0.0005402947  0.0005402947
## Soft        0.0007415283  0.0005402947 -0.0012818231
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##       Hard Semi-rigid       Soft 
##  0.3333333  0.3333333  0.3333333
## Done.
so3<-make.simmap(treenew,x,model="ARD",nsim=nsim[3])
## make.simmap is sampling character histories conditioned on
## the transition matrix
## 
## Q =
##                     Hard    Semi-rigid          Soft
## Hard       -0.0001997071  0.0000000000  0.0001997071
## Semi-rigid  0.0023041591 -0.0023041591  0.0000000000
## Soft        0.0021232656  0.0005497209 -0.0026729866
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##       Hard Semi-rigid       Soft 
##  0.3333333  0.3333333  0.3333333
## Done.
treenews<-c(so1, so2, so3)

objs<-describe.simmap(treenews)

plot(objs,type="fan",fsize=0.01,lwd=1,ftype="i",
     colors=cols,ylim=c(-2,Ntip(treenew)),part=0.97)
add.simmap.legend(colors=cols,prompt=FALSE,x=-300,y=280)

Same result as with semi-rigid coding, unsurprisingly.

  • Coded as hard
x[c(21,22)]<-"Hard"

logL<-sapply(c("ER","SYM","ARD"),
    function(model,treenew,x) make.simmap(treenew,x,model)$logL,
    tree=treenew,x=x)
## make.simmap is sampling character histories conditioned on
## the transition matrix
## 
## Q =
##                     Hard    Semi-rigid          Soft
## Hard       -0.0007710002  0.0003855001  0.0003855001
## Semi-rigid  0.0003855001 -0.0007710002  0.0003855001
## Soft        0.0003855001  0.0003855001 -0.0007710002
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##       Hard Semi-rigid       Soft 
##  0.3333333  0.3333333  0.3333333
## Done.
## make.simmap is sampling character histories conditioned on
## the transition matrix
## 
## Q =
##                     Hard   Semi-rigid          Soft
## Hard       -0.0006317588  0.000000000  0.0006317588
## Semi-rigid  0.0000000000 -0.000568647  0.0005686470
## Soft        0.0006317588  0.000568647 -0.0012004058
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##       Hard Semi-rigid       Soft 
##  0.3333333  0.3333333  0.3333333
## Done.
## make.simmap is sampling character histories conditioned on
## the transition matrix
## 
## Q =
##                     Hard   Semi-rigid          Soft
## Hard       -0.0004499743  0.000000000  0.0004499743
## Semi-rigid  0.0034859171 -0.007000114  0.0035141964
## Soft        0.0003928020  0.000000000 -0.0003928020
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##       Hard Semi-rigid       Soft 
##  0.3333333  0.3333333  0.3333333
## Done.
AIC<-mapply(aic,logL,c(1,3,6))
AIC.W<-aic.w(AIC)

nsim<-1000
Nsim<-round(nsim*AIC.W)
d<-if(sum(Nsim)>nsim) -1 else 1
nsim<-Nsim+d*sample(c(rep(1,abs(nsim-sum(Nsim))),
    rep(0,length(Nsim)-abs(nsim-sum(Nsim)))))
nsim
##  ER SYM ARD 
## 347 497 156
h1<-make.simmap(treenew,x,model="ER",nsim=nsim[1])
## make.simmap is sampling character histories conditioned on
## the transition matrix
## 
## Q =
##                     Hard    Semi-rigid          Soft
## Hard       -0.0007710002  0.0003855001  0.0003855001
## Semi-rigid  0.0003855001 -0.0007710002  0.0003855001
## Soft        0.0003855001  0.0003855001 -0.0007710002
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##       Hard Semi-rigid       Soft 
##  0.3333333  0.3333333  0.3333333
## Done.
h2<-make.simmap(treenew,x,model="SYM",nsim=nsim[2])
## make.simmap is sampling character histories conditioned on
## the transition matrix
## 
## Q =
##                     Hard   Semi-rigid          Soft
## Hard       -0.0006317588  0.000000000  0.0006317588
## Semi-rigid  0.0000000000 -0.000568647  0.0005686470
## Soft        0.0006317588  0.000568647 -0.0012004058
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##       Hard Semi-rigid       Soft 
##  0.3333333  0.3333333  0.3333333
## Done.
h3<-make.simmap(treenew,x,model="ARD",nsim=nsim[3])
## make.simmap is sampling character histories conditioned on
## the transition matrix
## 
## Q =
##                     Hard   Semi-rigid          Soft
## Hard       -0.0004499743  0.000000000  0.0004499743
## Semi-rigid  0.0034859171 -0.007000114  0.0035141964
## Soft        0.0003928020  0.000000000 -0.0003928020
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
##       Hard Semi-rigid       Soft 
##  0.3333333  0.3333333  0.3333333
## Done.
treenewh<-c(h1, h2, h3)

objh<-describe.simmap(treenewh)

plot(objh,type="fan",fsize=0.01,lwd=1,ftype="i",
     colors=cols,ylim=c(-2,Ntip(treenew)),part=0.97)
add.simmap.legend(colors=cols,prompt=FALSE,x=-300,y=280)

Ancestral state for dinosaurs and archosaurs, and all major less inclusive clades except pterosaurs, is hard-shelled, as in the Legendre et al. (2020) coding.

To check how strongly are these results influenced by branch length information, we can replicate these analyses using maximum parsimony, which does not consider branch length information, using castor.

  • Coding from this study
studyN<-thisstudy
studyN[studyN=="Soft"]<-1; studyN[studyN=="Semi-rigid"]<-2; studyN[studyN=="Hard"]<-3
studyN<-as.numeric(studyN); names(studyN)<-names(thisstudy)

MPS<-asr_max_parsimony(tree, studyN, Nstates=3)
plotTree(tree,type="fan",fsize=0.01,lwd=1,ftype="i",
         colors=cols,ylim=c(-2,Ntip(treenew)),part=0.97)
nodelabels(node=1:tree$Nnode+Ntip(tree),pie=MPS$ancestral_likelihoods,piecol=cols,cex=0.5)
tiplabels(pie=to.matrix(studyN,sort(unique(studyN))),piecol=cols,cex=0.3)
add.simmap.legend(colors=cols,prompt=FALSE,x=-300,y=280,fsize=0.8)

# ancestral states
MPS2<-MPS$ancestral_likelihoods
rownames(MPS2)<-c(209:374); colnames(MPS2)<-c("soft","semi-rigid","hard")
  • Coding from Norell et al. (2020)
norellN<-norell
norellN[norellN=="Soft"]<-1; norellN[norellN=="Semi-rigid"]<-2; norellN[norellN=="Hard"]<-3
norellN<-as.numeric(norellN); names(norellN)<-names(norell)

MPN<-asr_max_parsimony(tree, norellN, Nstates=3)
plotTree(tree,type="fan",fsize=0.01,lwd=1,ftype="i",
         colors=cols,ylim=c(-2,Ntip(treenew)),part=0.97)
nodelabels(node=1:tree$Nnode+Ntip(tree),pie=MPN$ancestral_likelihoods,piecol=cols,cex=0.5)
tiplabels(pie=to.matrix(norellN,sort(unique(norellN))),piecol=cols,cex=0.3)
add.simmap.legend(colors=cols,prompt=FALSE,x=-300,y=180,fsize=0.8)

# ancestral states
MPN2<-MPN$ancestral_likelihoods
rownames(MPN2)<-c(209:374); colnames(MPN2)<-c("soft","semi-rigid","hard")
  • Coding from Legendre et al. (2020)
legendreN<-legendre
legendreN[legendreN=="Soft"]<-1; legendreN[legendreN=="Semi-rigid"]<-2; legendreN[legendreN=="Hard"]<-3
legendreN<-as.numeric(legendreN); names(legendreN)<-names(legendre)

MPL<-asr_max_parsimony(tree, legendreN, Nstates=3)
plotTree(tree,type="fan",fsize=0.01,lwd=1,ftype="i",
         colors=cols,ylim=c(-2,Ntip(treenew)),part=0.97)
nodelabels(node=1:tree$Nnode+Ntip(tree),pie=MPL$ancestral_likelihoods,piecol=cols,cex=0.5)
tiplabels(pie=to.matrix(legendreN,sort(unique(legendreN))),piecol=cols,cex=0.3)
add.simmap.legend(colors=cols,prompt=FALSE,x=-300,y=280,fsize=0.8)

# ancestral states
MPL2<-MPL$ancestral_likelihoods
rownames(MPL2)<-c(209:374); colnames(MPL2)<-c("soft","semi-rigid","hard")
  • With the reduced tree, changing only the two sauropodomorphs
xN<-x
xN[xN=="Soft"]<-1; xN[xN=="Semi-rigid"]<-2; xN[xN=="Hard"]<-3
xN<-as.numeric(xN); names(xN)<-names(x)

# Coded as soft
xN[c(21,22)]<-1
MPx<-asr_max_parsimony(treenew, xN, Nstates=3)
plotTree(treenew,fsize=0.4,lwd=1,ftype="i",colors=cols)
nodelabels(node=1:treenew$Nnode+Ntip(treenew),pie=MPx$ancestral_likelihoods,piecol=cols,cex=0.5)
tiplabels(pie=to.matrix(xN,sort(unique(xN))),piecol=cols,cex=0.3)
add.simmap.legend(colors=cols,prompt=FALSE,x=25,y=180,fsize=0.8)

MPx2<-MPx$ancestral_likelihoods
rownames(MPx2)<-c(202:360); colnames(MPx2)<-c("soft","semi-rigid","hard")

# Coded as semi-rigid
xN[c(21,22)]<-2
MPx<-asr_max_parsimony(treenew, xN, Nstates=3)
plotTree(treenew,fsize=0.4,lwd=1,ftype="i",colors=cols)
nodelabels(node=1:treenew$Nnode+Ntip(treenew),pie=MPx$ancestral_likelihoods,piecol=cols,cex=0.5)
tiplabels(pie=to.matrix(xN,sort(unique(xN))),piecol=cols,cex=0.3)
add.simmap.legend(colors=cols,prompt=FALSE,x=25,y=180,fsize=0.8)

MPx2<-MPx$ancestral_likelihoods
rownames(MPx2)<-c(202:360); colnames(MPx2)<-c("soft","semi-rigid","hard")

# Coded as hard
xN[c(21,22)]<-3
MPx<-asr_max_parsimony(treenew, xN, Nstates=3)
plotTree(treenew,fsize=0.4,lwd=1,ftype="i",colors=cols)
nodelabels(node=1:treenew$Nnode+Ntip(treenew),pie=MPx$ancestral_likelihoods,piecol=cols,cex=0.5)
tiplabels(pie=to.matrix(xN,sort(unique(xN))),piecol=cols,cex=0.3)
add.simmap.legend(colors=cols,prompt=FALSE,x=25,y=180,fsize=0.8)

MPx2<-MPx$ancestral_likelihoods
rownames(MPx2)<-c(202:360); colnames(MPx2)<-c("soft","semi-rigid","hard")

Ancestral state reconstruction – continuous character (eggshell thickness)

  • Data and tree
treedi<-multi2di(tree, random=FALSE)

# To visualize node numbers on the tree
ggtree(treedi) + geom_text2(aes(subset=!isTip, label=node), hjust=-.3) + geom_tiplab()

AET<-log1p(datanew$Eggshell_thickness); names(AET)<-rownames(datanew)
RET<-log1p(datanew$Eggshell_thickness/datanew$Egg_mass); names(RET)<-rownames(datanew)
  • Phylogenetic signal (contMap assumes a Brownian motion model)
phylosig(treedi, AET, method="lambda", test=TRUE)
## 
## Phylogenetic signal lambda : 0.999934 
## logL(lambda) : -316.25 
## LR(lambda=0) : 255.021 
## P-value (based on LR test) : 2.08867e-57
phylosig(treedi, RET, method="lambda", test=TRUE)
## 
## Phylogenetic signal lambda : 0.999934 
## logL(lambda) : -231.389 
## LR(lambda=0) : 206.653 
## P-value (based on LR test) : 7.38001e-47

Very strong signal for both absolute and relative eggshell thickness.

  • ASR for absolute eggshell thickness
phyloplotAET<-contMap(treedi, AET, plot=F)
plot(setMap(phyloplotAET,colors=rev(brewer.pal(10, "Spectral"))),type="fan",
     fsize=0.01,lwd=4,ylim=c(-2,Ntip(treenew)),part=0.8,legend=FALSE)
add.color.bar(leg=100,cols=rev(brewer.pal(10,"Spectral")),
              prompt=FALSE,x=-300,y=280)

# Ancestral states
expm1(fastAnc(treedi, AET))
## Ancestral character estimates using fastAnc:
##         209         210         211         212         213         214 
##   30.237176   30.237176   27.637279   20.146957   37.147079   57.434822 
##         215         216         217         218         219         220 
##   59.490746   38.962770   37.947764   62.376276   63.738529   60.890954 
##         221         222         223         224         225         226 
##   18.797940   17.535632   17.058187   21.413045   14.494621   17.213248 
##         227         228         229         230         231         232 
##   16.971098   16.477399   14.535149   16.143759    8.220722    8.894198 
##         233         234         235         236         237         238 
##   17.611084   17.663373   11.373368   11.394027   11.094212    7.685149 
##         239         240         241         242         243         244 
##   18.197125   27.804160   11.643952    6.395467    5.262258   11.872530 
##         245         246         247         248         249         250 
##   10.876488   12.263387   12.123808   14.013434   13.895900   14.226852 
##         251         252         253         254         255         256 
##   16.558085   31.172635   40.914364  133.738678  544.568527  573.610635 
##         257         258         259         260         261         262 
##  610.795533  607.951060  498.450624  137.074853  142.357447  143.074080 
##         263         264         265         266         267         268 
##  139.614328  143.141126   61.052662  155.004338  178.337602  198.598257 
##         269         270         271         272         273         274 
##  196.829154  167.286930  145.639082  217.263595  259.777176  324.061138 
##         275         276         277         278         279         280 
##   38.519274  253.555625  279.270796  414.658851  624.034336  367.028433 
##         281         282         283         284         285         286 
##  460.030533  450.070207  468.619867  478.051897   33.579247    2.962142 
##         287         288         289         290         291         292 
##    0.014681    4.087266   11.285567    3.380069   41.916492  152.227600 
##         293         294         295         296         297         298 
##  152.227600  152.227600 1426.916211 1426.916211  152.227600 1818.852274 
##         299         300         301         302         303         304 
## 1818.852274 1455.301611 1501.790647   31.092007   14.703935   72.646286 
##         305         306         307         308         309         310 
##    1.813370 1725.926348 1725.926348 1725.926348 1725.926348 1725.926348 
##         311         312         313         314         315         316 
## 1725.926348 1725.926348 1725.926348 1725.926348 1725.926348 1725.926348 
##         317         318         319         320         321         322 
## 1725.926348 1725.926348 1725.926348 1725.926348 1725.926348 1725.926348 
##         323         324         325         326         327         328 
## 1725.926348  321.949973  479.318439  856.292119  856.292119  856.292119 
##         329         330         331         332         333         334 
##  533.581298  535.746398  535.746398  557.800038 1286.491850  567.795787 
##         335         336         337         338         339         340 
##  606.755787  606.755787  606.755787  606.755787 1391.768701 1101.675921 
##         341         342         343         344         345         346 
##  540.855621  540.855621  540.855621  656.650226  418.648969  418.648969 
##         347         348         349         350         351         352 
##  418.648969  400.036092  400.036092  400.036092  400.036092  383.759407 
##         353         354         355         356         357         358 
##  393.341060  482.609908  482.609908  482.609908  482.609908  482.609908 
##         359         360         361         362         363         364 
##  487.383928  482.609908  672.440036  688.378914  743.268097  749.865095 
##         365         366         367         368         369         370 
##  642.582251 1304.089170  319.366201  778.769338  908.526782  820.791309 
##         371         372         373         374         375         376 
##  423.102753  255.042210  300.656066  343.613100  298.832008  290.376178 
##         377         378         379         380         381         382 
##  222.124964  217.737359  177.855640  166.738135  161.028661  229.680422 
##         383         384         385         386         387         388 
##  234.234604  236.215710  233.547118  232.188904  231.149423  196.464412 
##         389         390         391         392         393         394 
##  216.848165  210.911249  208.760967  210.718962  206.849992  201.750368 
##         395         396         397         398         399         400 
##  251.215833  289.694249  316.988394  253.496044  257.536715  243.363765 
##         401         402         403         404         405         406 
##  151.833016  155.871496  164.021661  144.782394  128.404579  106.959834 
##         407         408         409         410         411         412 
##  100.296809  118.839823  109.667780   73.949894  103.667666   99.285850 
##         413         414         415 
##   98.490499  103.797577  101.497881

Ancestral eggshell thickness seems to be intermediate for all major nodes, including archosaurs and dinosaurs. However, many extant lepidosaurs are also recovered as intermediate – likely due to the extremely low values in pterosaurs, reported as lacking a calcareous layer entirely, which shift any thicker eggshell towards the middle of the spectrum.

Furthermore, the pattern strongly follows that of body mass, as expected (e.g. Stein et al., 2019; Legendre and Clarke, 2021).

  • ASR for relative eggshell thickness
phyloplotRET<-contMap(treedi, RET, plot=F)
plot(setMap(phyloplotRET,colors=rev(brewer.pal(10, "Spectral"))),type="fan",
     fsize=0.01,lwd=4,ylim=c(-2,Ntip(treenew)),part=0.8,legend=FALSE)
add.color.bar(leg=100,cols=rev(brewer.pal(10,"Spectral")),
              prompt=FALSE,x=-300,y=280)

# Ancestral states
expm1(fastAnc(treedi, RET))
## Ancestral character estimates using fastAnc:
##       209       210       211       212       213       214       215       216 
##  3.463807  3.463807  4.072394  7.776137 20.371286 37.945868 39.830943 25.691073 
##       217       218       219       220       221       222       223       224 
## 39.657850 42.685185 44.863801 53.013755  7.517838  7.370333  7.179543  6.555069 
##       225       226       227       228       229       230       231       232 
##  9.463869  7.389222  7.661106  9.194875 11.651921 13.025134 12.735303 13.646880 
##       233       234       235       236       237       238       239       240 
## 16.532983 20.147667 21.875358 24.370764 23.821319 14.769232  1.647382  1.003291 
##       241       242       243       244       245       246       247       248 
##  1.610748  0.936602  0.860325  1.828704  1.596699  2.012611  3.474841  1.683005 
##       249       250       251       252       253       254       255       256 
##  1.180153  1.164833  2.085560  3.133937  2.664801  6.649956 10.831692 10.811365 
##       257       258       259       260       261       262       263       264 
## 10.536241 10.641737  9.256815  6.921701  8.109834 10.716757 10.797700 12.176394 
##       265       266       267       268       269       270       271       272 
##  2.043793 11.104335 10.749829 10.936284 10.273487  9.993487  8.767791 10.221733 
##       273       274       275       276       277       278       279       280 
##  7.286859  7.960865  2.192523  4.759820  4.819961  5.480688  6.780228  4.590287 
##       281       282       283       284       285       286       287       288 
##  5.810148  6.065964  5.670210  5.486058  1.755410  0.464691  0.004048  0.554308 
##       289       290       291       292       293       294       295       296 
##  1.153786  0.409454  1.322077  1.359789  1.359789  1.359789  5.337305  5.337305 
##       297       298       299       300       301       302       303       304 
##  1.359789  0.659302  0.659302  1.233747  1.279996  0.923542  0.437393  0.389384 
##       305       306       307       308       309       310       311       312 
##  0.131437  0.730699  0.730699  0.730699  0.730699  0.730699  0.730699  0.730699 
##       313       314       315       316       317       318       319       320 
##  0.730699  0.730699  0.730699  0.730699  0.730699  0.730699  0.730699  0.730699 
##       321       322       323       324       325       326       327       328 
##  0.730699  0.730699  0.730699  2.798346  4.269776  4.207026  4.207026  4.207026 
##       329       330       331       332       333       334       335       336 
##  5.096371  5.317186  5.317186  5.249324  2.804855  5.265575  4.919952  4.919952 
##       337       338       339       340       341       342       343       344 
##  4.919952  4.919952  2.041476  3.280823  5.648622  5.648622  5.648622  4.639335 
##       345       346       347       348       349       350       351       352 
##  8.397269  8.397269  8.397269  8.721975  8.721975  8.721975  8.721975  8.944506 
##       353       354       355       356       357       358       359       360 
##  8.806187  8.062933  8.062933  8.062933  8.062933  8.062933  8.107327  8.062933 
##       361       362       363       364       365       366       367       368 
##  3.124942  3.017798  2.342437  2.259183  2.383122  0.661661  4.634284  2.112206 
##       369       370       371       372       373       374       375       376 
##  1.859723  1.878858  1.176578 10.834138  9.311660  7.858083  9.853769 10.770714 
##       377       378       379       380       381       382       383       384 
## 11.860656 11.582529 13.304265 13.852859 14.336956 10.818602 10.468547 10.019546 
##       385       386       387       388       389       390       391       392 
##  9.677966  9.355086  9.031513 13.464317 11.100694 10.828751 10.901636 10.572284 
##       393       394       395       396       397       398       399       400 
## 10.469217 10.385413  9.574631  8.311346  7.455384  9.415512  8.811797  9.272837 
##       401       402       403       404       405       406       407       408 
## 19.516988 20.275505 19.452318 20.563944 23.256399 28.370771 30.397447 25.081937 
##       409       410       411       412       413       414       415 
## 27.482967 41.577914 29.373696 30.524621 32.153419 31.108670 34.366359

Much lower values for archosaurs and dinosaurs – eggshell thickness seems to have become thinner for a given egg mass at the base of Ornithodira, and later increased in theropods. However, this pattern is dependent on a very small sample of pterosaurs, ornithischians, and sauropods; it is likely that the addition of new specimens attributed to any of these clades would considerably change that pattern.

We can see a strong increase in geckos and in eufalconimorphs – the latter having already been identified in Legendre and Clarke (2021).

References

  • Legendre, L.J., Clarke, J.A., 2021. Shifts in eggshell thickness are related to changes in locomotor ecology in dinosaurs. Evolution 75, 1415–1430.
  • Legendre, L.J., Rubilar-Rogers, D., Musser, G.M., Davis, S.N., Otero, R.A., Vargas, A.O., Clarke, J.A., 2020. A giant soft-shelled egg from the Late Cretaceous of Antarctica. Nature 583, 411–414.
  • Norell, M.A., Wiemann, J., Fabbri, M., Yu, C., Marsicano, C.A., Moore-Nall, A., Varricchio, D.J., Pol, D., Zelenitsky, D.K., 2020. The first dinosaur egg was soft. Nature 583, 406–410.
  • Stein, K., Prondvai, E., Huang, T., Baele, J.-M., Sander, P.M., Reisz, R., 2019. Structure and evolutionary implications of the earliest (Sinemurian, Early Jurassic) dinosaur eggs and eggshells. Scientific Reports 9, 4424.